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PREFACE 

This  report  provides  an  overview  of  some  advanced  digital  signal  processing  techniques  for 
onc'  and  two-dimensional  data  that  may  be  useful  for  pattern  recognition,  real-time  image  process¬ 
ing,  and  on-line  test  article  health  monitoring  functions.  Topics  include  a  review  of  classical  Fourier 
analysis  and  its  application  to  transient  data,  the  Short  Time  Fourier  Transform,  multi-resolution 
analysis,  wavelet  analysis,  and  applications  of  wavelet  analysis  to  some  typical  AEDC-typ)e  signals. 

The  work  reported  herein  was  conducted  at  Arnold  Engineering  Development  Center 
(AEDC),  Air  Force  Materiel  Command  (AFMC),  under  Program  Element  65807F  at  the  request  of 
AEDC/DOT,  Arnold  Air  Force  Base,  TN.  The  AEDC/IX)T  Project  Manager  was  James  D. 
Mitchell.  Management  for  this  project  was  performed  by  Sverdrup  Technology,  Inc.,  AEDC  Group, 
technical  services  contractor  of  the  propulsion  test  facilities,  AEDC,  AFMC,  Arnold  Air  Force 
Base,  TN,  under  Air  Force  Project  No.  DC97EW  (Job  0088).  The  Sverdrup  Project  Manager  was 
T.  F.  Tibbals.  This  work  was  performed  in  collaboration  with  the  University  of  Tennessee  Space 
Institute,  Dr.  K.  R.  Kimble,  principal  investigator,  under  contract  A94W-07  with  Sverdrup 
Technology,  Inc.  The  manuscript  was  submitted  for  publication  on  October  31, 1994. 
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1.0  INTRODUCTION 

Various  dynamic  data  analysis  tools  for  one^imensional  time-domain  signals  are  employed 
to  determine  the  Irequency  content  of  a  signal  for  mechanical  analysis.  When  it  is  tied  to  the 
fundamental  frequency  of  the  various  components  comprising  the  machinery  being  evaluated,  this 
information  gives  an  indication  of  the  state  or  health  of  the  machine.  Current  techniques  for 
evaluating  dynamic  data  for  potential  mechanical  problems  are  primarily  centered  around  the  Fast 
Fourier  Transform  (FFT)  and  the  Short  Time  Fourier  Transform  (STFT).  However,  the  use  of 
Fourier  analysis  for  frequency  component  extraction  is  restricted  to  bandlimited  stationary  signals. 
Because  of  this  restriction,  small  transients  may  not  be  detected  due  to  a  smoothing  effect  of  the 
FFT,  or  the  FFT  spectrum  may  be  smeared  due  to  frequency  ramping  and  abrupt  incidents  or 
discontinuities  in  the  signal.  Various  techniques  have  been  employed  to  overcome  the  limitations 
of  the  FFT  for  non-stationary  data.  These  techniques  include  windowed  Fourier  Transform  (Gabor 
or  Short  Time  Fourier  Transform),  synchronous  sampling  to  remove  RPM-ramp  effects,  Wigner- 
Ville  analysis,  and,  more  recently,  wavelet  analysis.  This  report  provides  an  overview  of  Fourier 
analysis  (including  the  STFT)  and  a  background  in  wavelet  theory  based  upon  the  analyzing 
function  basis  approach  so  that  wavelet  theory  may  be  contrasted  against  Fourier  analysis.  A  model 
signal  with  stationary  and  transient  characteristics  is  developed  to  permit  comparisons  of  various 
analysis  techniques  based  upon  a  known  analytic  signal  which  resembles  a  real  vibration  signal. 
Some  applications  of  wavelets  to  other  transient  signals  are  also  provided. 

2.0  FOURIER  TRANSFORM 


Consider  the  representation  of  a  finite-power  signal,  40,  defined  on  the  interval  (tQ,  %  +  T)  in 

terms  of  a  set  of  preselected  time  functions  d>|(0 . K  is  convenient  to  choose  these 

functions  with  properties  analogous  to  the  orthogonal  unit  vectors  of  Cartesian  space,  i.e., 
orthogonal  such  that 

/  $„(<)$;(t)dt  =  0  Vm  #  n 

Ju 

and  normalized  such  that 


i 


=  I 


where  <I>*  (/)  is  the  complex  conjugate  of  (r).  Therefore,  these  functions  are  orthonormal 
functions  on  the  interval  (/q,  Iq  -i-  7),  where  the  orthonormality  condition  was  chosen  since  the 
signal,  being  periodic,  is  bounded  in  power  though  unbounded  in  energy. 
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Now  consider  an  approximation  x  (/)  to  x( t)  based  upon  a  series  expansion  of  the  form 

N 

n=l 


where  is  a  given  finite  positive  integer  and  where  the  constants  C„  are  to  be  chosen  such  that  x  (/) 
represents  x( t)  as  closely  as  possible.  A  criterion  of  measure  for  this  closeness  of  approximation 
must  be  chosen  based  upon  the  signal's  characteristics.  In  this  case,  we  choose  to  measure  the  error 
in  a  square-integral  sense  since  our  concern  is  with  power  signals.  We  define  this  error  as 

tAT  =  /  |a;(0  -  . 

■'to 


and  we  wish  to  find  the  C„  coefficients  to  minimize  this  error.  Substituting  for  x  (r)  yields 


Jrlo+T 
to 


N 


n=l 

^  J  fto+T  yio+T 

Ktrdi-J2\C'„  x(t)'l^-„{t)di  +  Cn 

■'‘k  n=l  L  ■^*0  dto 

Adding  and  subtracting  ^  (f„j 


n=l 


I  +r 

r 

'b 

N  1  /-Co-hT 


dt  yields 


fio+T  N  ,5+t  JV  <,+!■ 

Jlo  „-J  ■'to  -'*0 


di  . 


Now  since  the  first  two  terms  above  are  independent  of  the  C„  coefficients  of  our  approxima¬ 
tion,  and  the  last  term  is  non-negative  and  adds  to  these  two  terms,  then  the  only  way  to  minimize 
the  error  is  to  choose  the  approximation  coefficients  C„  such  that 

yto+^ 

C„=  x[i)K{l)dl  n=l,2,...,N 

•fta 

Note  that  this  is  an  inner  product  such  that  C„  =  <x,  0„  >.  For  special  choices  of  the  set  of  functions 

(t>l(r) . referred  to  as  complete  sets  in  the  Hilbert  function  space  I?  of  all  square-integrable 

functions,  it  will  be  true  that 


lim  tN.min  =  iim 

A— ^00  jV— oo 


n--M 


=  0 
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for  any  signal  that  is  square-integrable,  Le., 


dt  <  oo  . 


In  the  case  of  zero  square-integral  error. 


and 


n=I 

■*  “''o  fi=:l 


This  is  the  generalized  Parseval  theorem  which  states  that  the  mean  power  of  the  periodic  signal  is 
the  sum  of  the  squares  of  the  approximating  coefficients.  The  exponential  functions  of  the  Fourier 
series  form  a  complete  orthogonal  set  in  I?  such  that 


^n{i) 


anuot 
^  ) 


n  =  0,±1,  ±2,  ...,oo 


are  axes  of  [?  and  Fourier  analysis  projects  jc(/)  onto  these  axes  just  as  in  the  case  of  breaking  a 
vector  into  its  orthogonal  components  is  to  project  the  vector  onto  the  orthogonal  axes  of  Cartesian 
space.  Now  (i>„=  llcfT  v/heie  T  is  the  fundamental  period,  that  is,  T  is  the  minimum  period  for  a 
complete  sine  wave  of  the  lowest  frequency  of  interest.  Thus,  partial  sums  of  exponential  Fourier 
series  minimize  the  square-integrable  error  between  the  series  and  the  signal  under  investigation 
(Ref.  1). 


Since  sine  waves  have  a  period  of  2t[  ,  a  function  x(r)  that  is  square  integrable  (i.e.,  piecewise 
continuous)  on  the  interval  (0, 2n)  i.e.. 


is  an  element  of  the  Hilbert  function  space  £.^(0,  271)  of  all  square  integrable  functions  defined  on 
(0,  2ti)  and  expansions  and  contractions  of  x(0  (i.e..  varying  the  number  of  cycles  on  the  interval 
(0,  2r))  form  a  basis  of  the  2ic)  function  space.  By  basis  is  meant  a  set  of  linearly  independent 
functions  which  can  be  combined  to  form  another  function  which  also  exists  in  the  same  function 
space.  Then  x(t)  has  a  Fourier  series  representation: 

OO 

x{i)=  2 

OO 
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where,  for  n  =  -«». . -2,  -1, 0, 1,2,.. .«», 


Cr  — 


Notice  that  the  number  of  basis  functions  is  infinite  for  this  representation. 
Per  Chui  (Ref.  2),  this  series  converges  in  l}  such  that 


Jim  „ 
Ar,.v-,oo  2-f 


r  /  WO-  E  =  0 

•'0  n=-M 


There  are  three  distinct  features  of  x{t)  in  the  Fourier  series  representation.  First,  x(0  is 
decomposed  into  mutually  orthogonal  components  since  the  inner  product 


2ir 


m  ^  n 
m  =  n 


Second,  the  bases  of  this  orthonormal  representation  are  all  generated  from  a  single  function  &)(/) 
=  e'®"'  by  integral  dilations,  i.e.,  w„(0  =  w(nt)  for  all  integers  n.  Finally,  the  Fourier  expansion 
coefficients  are  determined  via  an  inner  product  =<x,  h'„>.  “It  is  orthogonality  that  allows  us 
to  find  each  term  separately  [no  cross-terms  means  no  cross-talk  or  interdependence  between 
frequencies],  and  it  is  completeness  that  allows  the  sines  and  cosines  to  reproduce  jr(r).”  (Ref.  3). 

However,  the  use  of  this  complex  exponential  basis  is  not  without  complications.  First,  the 
approximation  using  Fourier  coefficients  is  an  infinite  summation  which  can  only  be  made  finite  by 
bandlimiting  the  signal.  Second,  the  basis  functions  themselves  have  infinite  extent  in  the  time 
domain,  and  therefore  cannot  readily  approximate  a  short-lived  event  in  the  signal.  Third,  the  signal 
must  be  periodic,  or  at  least  the  results  of  Fourier  analysis  assume  periodicity,  which  may  not  be 
the  case  for  typical  rea/  signals.  Finally,  the  signal  should  meet  the  Dirichlet  conditions  to  ensure 
the  existence  of  the  Fourier  series  approximation  *  (i.e..  the  series  converges  to  the  continuous  signal 
x(t}}.  The  Dirichlet  conditions  are: 


1.  Note  that  if  these  conditions  hold  as  T  — »  <»,i.  e.,  the  signal  is  aperiodic^  then  the  signal’s  Fourier 
Integra]  Transfonn  exists  and  the  signal  can  be  approximated  by  a  sum  of  pure  sinusoids  which  vary 
continuously  in  frequency.  In  this  case  the  spectrum  is  continuous  instead  of  discrete  as  in  the  Fourier  Series 
representation. 
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•  If  the  function  has  discontinuities,  their  number  must  be  finite  in  any  period. 

•  The  function  must  contain  a  finite  number  of  maxima  and  minima  during  any  period. 

•  The  function  must  be  absolutely  integrable  in  any  period,  i.e., 

f  |x(<)|(ft  <  CX) 

Jo 

where  jr(0  is  a  continuous  (V(  at  least  piecewise)  describing  function  which  approximates  the 
actual  (real)  signal.  Note  that  the  Dirichlet  conditions  impose  stationarity  onto  the  signal.  In 
addition,  note  that  there  is  no  time  domain  translation  parameter,  as  will  be  seen  later  for  wavelets, 
since  it  makes  no  sense  to  translate  an  infinite  extent  function  —  this  is  handled  by  the  concept  of 
phase.  However,  this  implies  no  time  domain  localization  is  possible,  and  consequently,  the  Fourier 
coefficients  are  independent  of  time. 

For  a  sampled  (digital)  signal  x,  (nA/)  as  an  approximation  to  xif)  (where  Af  =  sample  interval 
and  T = NM  analysis  time  block),  we  can  then  approximate  as 

j  N-\ 

Cj,  =f  -  ^  (x,(nAt)e-'*"«’‘^‘)Il 

n=0  ^ 

j  N-l 

n=0 

W-I 

n=a 

where  t  — ^  nAf  and  A/  — >  T/N ,  and  N  is  the  number  of  samples  (measurement  values)  of  the  signal 
within  a  given  localized  lime  interval  of  the  sampling  window.  This  defines  the  Eliscrete  Fourier 
Transform  (DFT)  for  a  bandlimited  signal.  The  Fast  Fourier  Transform  (FFT)  is  merely  a  fast 
algorithm  for  performing  the  DFT. 

3.0  SHORT  TIME  FOURIER  TRANSFORM  (GABOR  TRANSFORM) 

To  better  approximate  transient,  short-lived  or  time-localized  phenomena,  we  mustier 
the  Fourier  analysis  beyond  that  of  just  the  simple  windows  caused  by  digitization  of  the  signal  or 
the  accumulation  of  N  samples  for  use  in  the  FFT.  One  approach  is  to  window  the  signal  in  time  to 
better  approximate  the  local  characteristics.  The  Short  Time  Fourier  Transform  (STFT,  or  Gabor 
Trarisform)  attempts  to  provide  this  localization.  The  STFT  of  signal  J((0  is  defined  as 
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STFTit.u,)  =  r  x((  -  <)5(C)e-‘"“^c/C 
7-00 


where  the  Gaussian  window  g(0,  {?  :  [0.  x  ]).  centered  about  t/2  is  defined  as 


<ry/2ir 


and  <T  (standard  deviation)  is  the  inflection  point  of  concavity  of  the  curve  such  that  large  o 
indicates  a  flat  curve  and  small  tr  indicates  a  peaked  curve.  For  comparison  purposes  with  wavelets 
later,  let  us  define  the  window  g(0  to  include  the  transformation  exponential  such  that 


<7(0  = 


g-inwt 


essentially  windows  the  transform  basis  functions.  Then  the  STFT  is  again  an  inner  product,  <x, 
>,  and  linearity  is  preserved.  Per  Chui  (Ref.  2),  the  time  window  width  is  Ag^^t)  =  -Jlc  while 
the  frequency  (spectral)  window  width  is  AggifU)  =  Jl  /o,  and  these  window  widths  arc  constant 
for  the  analysis  independent  of  frequency,  as  shown  in  the  STFT  calculation  on  the  next  page.  Note 
that  the  product  of  these  window  widths  is  a  constant  in  agreement  with  the  Heisenberg  Uncertainty 
Principle. 

I 

A 


-k - k - k — k — i - k — i i - A — 1- 


t 


STFT  CALCULATION  GRID 
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The  STFT  can  be  discretized  following  the  same  approach  as  for  the  Discrete  Fourier 
Transform.  For  a  sampled  (digital)  signal  ar,(nAr)  as  an  approximation  to  x{t)  (where  At  =  sample 
interval  and  T=  AfAt  =  analysis  time  block),  we  can  then  approximate  for  ^  =  mAt,  m<N,  as 


'V-l  „ 

S  ^  (r.(mAt  -  nA0y(mA0e-“"‘'’"'^')^ 

y  AT-l 

=  —  ^  r,(mAt  —  nAf)3(mAt)c~**‘*'®i^^ 

n=0 

fi=0 

where  again  t  — >«Ar  and  A/  — >  T/N ,  is  the  number  of  samples  (measurement  values)  of  the  signal 
within  a  given  localized  time  interval  of  the  sampling  window,  and  M  is  the  number  of  samples 
within  the  Gaussian  window  e-{m-M/2j^/2a^  s’ppp  actually  has 

two  windows  applied  in  the  time  domain.  The  first  window  is  the  blocking  of  the  data  itself,  just  as 
in  the  normal  DFT.  The  second  window  is  the  Gaussian  window  which  localizes  the  signal  in  time 
by  making  the  signal  look  as  if  it  has  compact  support.  However,  one  must  remember  that  this 
compact  support  is  assumed  to  be  periodic.  In  addition,  the  window  must  be  sufficiently  narrow 
such  that  the  signal  is  stationary  within  the  window  so  that  the  Fourier  Transform  is  applicable. 
Obviously,  the  STFT  will  work  for  a  narrowband  signal  such  that  an  appropriate  window  can  be 
satisfactorily  narrowed  without  adversely  affecting  the  lower  hequencies. 

3.1  APPLICATION  TO  MODEL  SIGNAL 


The  use  of  the  STFT  for  gear  fault  detection  via  vibrational  analysis  is  described  in  Ref.  4. 
Gear  fault  detection  is  an  inherently  transient  analysis  since  a  problem  with  a  gear  tooth  can  only 
be  perceived  during  the  time  the  tooth  is  in  mesh  with  other  components.  The  STFT  time-frequency 
output  for  such  an  analysis,  however,  requires  some  interpretation.  Following  the  example  of  Ref. 
4,  consider  a  model  signal  which  resembles  a  gear  fault  and  is  composed  of  two  sinusoids  with 
stationary  properties,  another  which  is  amplitude  modulated,  and  two  Gaussian-shaped  impulses  as 
described  by 


*(ra)  =  Ai  sin 


2ir/in 

N 


+  (  1  -H  j42  cos 


N  . 


A3  «m 


2ir/Dn 

~ir~ 


+  Ate 


for  0<n<Ar-l. 
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The  parameters  used  were:  N=  1024, /o=  10  Hz,  /i  =  80  Hz, /2=  128  Hz,  Aj  =  1.5,  A2  =  0.75,  A3 
=  1 .0,  A4  =  3.0,  A5  =  3.0,  Bi  =  64,  ^2  =  64,  Hi  =  320,  and  «2  =  576  and  the  sample  rate  at  2.56F^ 
=  327.68  Hz  (where  is  the  maximum  static  frequency  component).  Gaussian  impulses  were 
used  to  simulate  a  real  mechanical  transient  event  rather  than  delta  functions.  The  Gaussian 
impulses  have  more  than  one  sample  describing  the  event  as  opposed  to  the  single  sample  of  a  delta 
function.  This  suggests  a  more  teal  energy  distribution  in  time  for  a  mechanical  system,  since  few 
mechanical  systems  can  instantaneously  respond  and  decay  as  implied  by  a  delta  function.  Figure 
1  shows  the  signal  and  its  spectrum.  Note  that  there  is  no  indication  of  the  Gaussian- shaped 
impulses,  and  the  amplitude  modulated  signal  component  at  128  Hz  has  split  into  two  components 
centered  around  128  Hz  separated  by  twice  the  modulation  frequency  of  =  10  Hz. 

Applying  the  STFT  to  this  signal  with  a  window  width  of  32  points  and  no  overlap  results  in 
the  spectrogram  of  Fig.  2.  The  stationary  sinusoids  appear  as  straight  horizontal  bands  in  the  time- 
frequency  plot,  continuous  for  the  whole  time  interval.  Note  that  the  frequency  resolution  is  coarse 
due  to  poor  frequency  localization  at  the  low  number  of  FFT  points,  while  the  time  resolution  ade¬ 
quately  shows  the  amplitude  modulation  at  /  =  128,  but  scalloping  has  smeared  the  true  amplitudes. 
The  Gaussian  impulses  have  lost  their  time  dependence,  and  only  contribute  to  an  amplitude  mod¬ 
ulation  effect  for  the  10- Hz  component,  since  the  impulses  span  more  than  one  FFT  window. 

Now  if  the  window  width  is  changed  to  128  points  to  provide  better  frequency  localization, 
then  the  spectrogram  of  Fig.  3  results.  Note  that  the  time  resolution  has  worsened  in  an  attempt  to 
better  localize  the  frequency  components.  In  this  case,  the  Gaussian  impulses  are  visible  since  the 
lO-Hz  peak  is  sharper  due  to  narrower  frequency  bin  spacing.  The  low-frequency  smearing  for  time 
slices  three  and  five  are  the  Gaussian  impulses  at  n  =  320  and  n  =  576,  respectively.  However,  note 
that  the  amplitude  modulation  at  128  Hz  is  no  longer  discernible,  and  that  it  is  impossible  to 
determine  whether  the  spikes  at  118  Hz  and  138  Hz  are  real  components  or  a  lO-Hz  amplitude 
modulation  of  a  128-Hz  component.  To  improve  time  resolution  with  the  same  size  FPTs,  the  signal 
was  analyzed  with  overlapping  windows.  Figures  4  and  5  illustrate  the  signal  and  its  spectrum  for 
SO-percent  overlap  with  32-  and  128-point  FPTs,  respectively.  Note  that  there  is  no  change  in  the 
50-percent  overlap  spectrogram  for  the  128 -point  FFTs  as  compared  to  the  zero  overlap 
spectrogram  (Fig.  3).  This  is  due  to  the  signal  characteristics  being  of  smaller  scale  than  what  even 
the  50-percent  overlapping  provides.  However,  the  32-point  FFT  spectrogram  with  50-percent 
overlap  was  able  to  discern  the  Gaussian  peaks  which  were  not  discernible  with  zero  overlap  (Fig. 
2).  This  improvement  in  information  content  via  overlapping  windows  is  at  the  heart  of  the  Malvar 
wavelet  approach,  and  will  be  discussed  further  in  Section  5.0. 

Again,  a  constant  window  width  was  utilized  for  all  frequencies  processed  by  the  STET, 
regardless  of  overlapping.  As  we  shall  see  with  wavelets,  the  ability  to  vary  the  window  width  with 
frequency  will  improve  the  analysis  (Ref.  2).  Also  note  that  by  changing  window  widths  in  a  series  of 
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such  analyses,  we  can  enhance  various  signal  characteristics  with  subsequent  trade-offs  in  time- 
fpequency  localizations  for  various  components.  This  series  of  analyses  constitutes  a  Multi-Resolution 
Analysis  (MRA)  as  discussed  in  the  next  section.  A  multi-resolution  analysis  will  allow  us  to  adjust 
our  windows  in  both  the  time  and  frequency  domains  to  optimize  time-riequency  localization. 

4.0  DAUBECHIES  WAVELET  TRANSFORM 
4.1  THE  MULTI-RESOLUTION  CONCEPT 

Multi-Resolution  Analysis  is  a  processing  technique  that  adapts  to  the  frequency  range  of 
interest  to  optimize  the  resolution  in  both  the  time  and  frequency  domains.  Generally,  while  main¬ 
taining  the  constraints  of  the  Heisenberg  Uncertainty  Principle,  this  analysis  adapts  to  the  frequency 
such  that  time  resolution  becomes  arbitrarily  good  at  high  frequencies  while  fluency  resolution 
becomes  arbitrarily  good  at  low  liequencies. 

Generally  and  mathematically,  a  multi-resolution  analysis  consists  of  an  increasing  sequence 
of  successive  approximation  spaces  Vj,  ranging  from  coarse  to  fine,  which  are  closed  subspaces 
within  and  satisfy  (following  the  notation  of  Daubechies  (Ref.  S) 

{0}  -i.  ■  •  •  c  Vj  C  Vi  C  Vb  C  VLi  C  VLj  C - >  X^(R)  =>  NESTED 


U  Vj  =  =>  CLOSURE 

j€2 

P  Vj  =  {0}  =>  ORTHOGONALITY 

}€Z 

ar(t)  €  V;-  <=►  ar(o<)  G  V;+i,  a  >  1  =}►  MULTI  -  RESOLUTION. 

If  Pj  is  an  orthogonal  projection  operator  onto  Vj,  then  the  closure  property  ensures  that  a  suit¬ 
able  approximation  to  xfr)  can  be  ultimately  obtained,  i.e., 

/irn  Pjx{i)  =  3:(t)  Vx  G  . 

The  multi-resolution  aspect  of  these  ladder  spaces  is  due  to  the  last  of  the  above  properties,  where 
all  spaces  are  a  scaled  version  of  each  other  such  that  ultimately  they  are  scaled  versions  of  a  central 
subspace  V^,  i.e., 

x{t)  e  Vj  <=>  x(a^t)  €  Vo,  a  >  1. 

Here  the  notion  of  scale  is  that  of  cartography  in  that,  at  a  given  scale,  the  signal  is  replaced  by  a 
best  approximation  that  can  be  drawn  at  that  scale.  By  moving  from  a  coarse  scale  to  a  fine  scale, 
one  zooms  in  on  the  details  of  the  signal,  and,  hence,  a  more  exact  representation  of  the  original 
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signal.  Obviously,  the  original  sampled  signal  is  all  the  information  that  exists  and  this  then 
becomes  the  finest  resolution  in  information  that  is  possible;  i.e.,  all  further  processing  can  only  de¬ 
compose  the  signal  to  those  bits  of  information  contained  in  the  signal  without  creating  any  new 
information. 

Multi-resolution  analysis  can  be  thought  of  as  an  approximation  sequence  based  upon  succes¬ 
sive  decimation  filters  where  the  finest  time  resolution  is  the  original  (typically  oversampled) 
signal.  Subsequent  coarser  time  domain  approximations  can  then  be  obtained  by  repeated 
application  of  decimate-by-n  filtering,  where  every  n  samples  are  removed  in  the  time  domain.  In 
the  frequency  domain,  this  results  in  peeling  away  the  higher  frequencies  with  each  application  of 
the  decimation  filter,  since  the  Nyquist  criterion  must  be  satisfied.  Now  to  increase  low-frequency 
resolution,  one  simply  accumulates  more  time-domain  samples  to  perform  a  larger  frequency 
domain  transform  such  that  the  same  N-point  FFT  results  in  NI2  frequency  bins  distributed  over  a 
much  smaller  frequency  bandwidth  due  to  the  filtering  and  decimation  process  [Dr.  B.W.  Bomar, 
UTSI].  Consequently,  time  domain  resolution  is  sacrificed  at  the  expense  of  frequency  domain 
resolution. 

Notice  that  one  can  also  look  at  multi-resolution  analysis  as  a  contraction  and  expansion  (or 
dilation)  of  the  support  of  the  projection/approximation  basis  functions.  This  is  due  to  the  fact  that 
at  high  frequency,  best  may  be  determined  to  be  an  A^-sampled  FFT;  therefore,  the  time  analysis 
window  is  Ntg,  where  tg  is  the  sample  time  interval.  This  means  that  the  support  of  the  analysis  is 
much  narrower  for  the  higher  frequencies  than  for  the  lower  frequencies  accumulated  over  nNtg, 
where  an  n-decimation  was  used  for  the  same  constant  q  =  Af/f, 

4,2  WAVELET  THEORY 

Wavelet  theory  unifies  various  signal  processing  techniques  developed  independently  to 
overcome  the  problems  of  Fourier  analysis.  For  example,  multi-resolution  signal  processing,  image 
compression  sub-band  coding,  and  wavelet  series  expansions  are  all  considered  a  single  theory 
mathematically  (Ref.  6).  Wavelet  theory  provides  a  good  technique  for  signal  approximation  using 
scaled  basis  functions,  of  which  Fourier  analysis  is  just  one  application  (where  the  support  is 
infinite,  non-decaying,  and  the  basis  functions  are  sine  functions),  and  allows  one  to  examine  the 
effectiveness  of  the  approximation  on  deeper  mathematical  levels  systematically.  This  paper  is 
primarily  concerned  with  the  application  of  wavelets  to  non-stationary  one-dimensional  signals, 
and  therefore  will  limit  its  discussion  of  wavelets  to  a  signal  processing  point  of  view.  An  additional 
interest  is  in  the  use  of  wavelets  to  systematically  improve  the  frequency  resolution  of  low- 
frequency  components  embedded  within  broadband  signals,  whether  stationary  or  non-stationary, 
without  unnecessarily  burdening  the  computation  of  high-frequency  spectral  components.  This 
then  will  allow  real-time  low-frequency  pattern  recognition  embedded  in  broadband  signals  such 
as  required  for  bearing  fault  diagnosis  based  upon  engine  case  accelerometer  data  (Ref  7), 
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The  Fourier  Transform  spanned  the  space  (0, 2Jt)  where  the  elements  x(fy  satisfy 


di  <oo  . 


However,  for  transient  signals,  there  may  be  no  periodicity  (at  least  for  the  transient  component 
even  if  embedded  in  a  periodic  signal);  therefore,  we  must  look  for  a  transform  which  span 
Z?(9l)  where  the  elements  x(/)  satisfy 


dt  <  oo  . 


As  stated  by  Chui  (Ref.  2,  p.  3), 


Clearly,  the  two  function  spaces  [for  each  of  these  transforms]  L^(0,  2n)  and 

(3t)  are  quite  different.  In  particular,  since  (the  local  average  values  of) 
every  function  in  L  (91)  must  “decay”  to  zero  at  ±oo ,  the  sinusoidal  (wave) 
functions  do  not  belong  to  Z?(95).  In  fact,  if  we  look  for  “waves”  that 
generate  L^(9l),  these  waves  should  decay  to  zero  at  ±<»;  and  for  all  practi¬ 
cal  purposes,  the  decay  should  be  very  fast.  That  is,  we  look  for  small  waves, 
or  “wavelets'*,  to  generate  I?  (91). 


The  rapid  decay  of  the  wavelet  provides  the  localization  necessary  to  adapt  to  local  transients  in  a 
non-stationary  waveform.  The  fact  that  we  look  for  a  wave  or  oscillating  function  at  all  is  to  some¬ 
how  preserve  the  concept  of  frequency  for  vibration/harmonic  analysis,  and  we  also  desire  local¬ 
ization  in  the  frequency  domain  as  well. 


Now  we  still  want  these  wavelets  to  be  generated  from  a  single  function  just  as  in  the  Fourier 
basis  functions  where 


u>n(t)  =  V)[nt)  =  V  integers  n 

to  ease  the  analysis.  In  addition,  we  still  want  the  functions  to  be  complete  and  orthogonal  in 
l}  (91)  so  we  can  be  assured  that  the  approximation  converges  to  x(j)  in  the  limit  and  that  we  can 
successively  find  each  coefficient  of  the  approximation. 

To  span  ]}  (91)  using  a  single  function  which  decays  rapidly,  i.e.,  is  compactly  supported  on 
the  real  line,  it  is  necessary  to  shift  the  function  along  91.  Consider  \|r  to  be  a  wavelet  basis  func¬ 
tion,  real,  and  of  compact  support  on  91;  then  for  t|/  to  span  all  of  91, 

^t{a:)  =  ^(a:  -  fr)  a:,6e». 
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This  is  the  translation  property  of  \|/ .  Now,  obviously  y  must  still  be  capable  of  describing  different 
frequencies,  since  our  intent  is  to  perform  a  frequency  analysis  of  the  mechanical  phenomena.  To 
accomplish  this  easily,  without  the  use  of  sinusoids,  we  must  consider  the  fundamental  definition 
of  frequency  as  the  inverse  of  the  period  between  zero  crossings  or  the  time  extent  of  the  compactly 
supported  analysis  functions.  With  this  in  mind,  we  consider  frequency  as  embedded  in  the  dilations 
of  the  basis  function  such  that 


^  -b)  I,  a,  t  e  3ft,  a  >  0 


Observe  that  our  development  parallels  that  for  Fourier  analysis,  and  seeks  a  basis  function 
that  is  obtained  from  a  single  function  by  dilations  (at)  to  adjust  to  varying  frequencies  and 
translations  (h/a)  to  span  the  real  line,  i.e.. 


«  ip(ax  —  6)  ft:  ^ 


i:,a,6  G  Sft,  a  >  0, 


and  is  different  from  the  Fourier  analysis  only  in  the  basis  functions  used,  which  are  not  periodic 
and  are  of  compact  support. 

Now  to  normalize  this  basis  in  consider 


II  'kiax  - 


Therefore 


II  ||2=a  5  |[^(j;)  II2 

so  that  for  be  of  unit  length, 

il>a,b{x)  =  a^il>{ax  —  b)  x,a,6e3ft,  a>0  . 

For  economy  of  the  analysis,  let's  restrict  ourselves  to  integral  dilations  and  translations  just  as  we 
restrict  ourselves  to  integral  frequencies  for  the  FFT.  LetZ  denote  the  set  of  integers;  then  for  j  €2, 
a-^aj-  t^,  >  1,  where  is  the  fixed  dilation  parameter.  Now  since  the  width  (support)  of  the 

analyzing  wavelet  will  be  proportional  to  it  follows  that  we  will  need  to  translate  in  steps 
proportional  to  the  width  of  the  wavelet  basis  function  to  cover  the  whole  real  line.  Consequently, 
let  ^0^1  ■>  ksZ,  where  b^ is  fixed.  Then 

-  kb^a-i) 

=  tij.)  leSft,  o>l,  60  >0,  j,k£Z. 
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Later,  will  become  the  sampling  rate;  however,  for  now  let  =  1,  where  it  is  asgnmp^  that  the 
shift  k  is  proportional  to  the  width  of  y . 

Since  our  interest  is  in  the  analysis  of  signals  evolving  in  time,  first  let  us  look  at  what  the  time 
description  of  the  wavelet  is.  Now  time  is  continuous  just  like  the  real  line.  However,  one  still  wants 
to  perform  integral  dilations  and  translations  utilizing  the  set  of  integers  Z  which  is  a  subset  of  91 . 
Therefore,  the  wavelets  describing  a  time  evolving  signal  instead  of  a  spatially  evolving  one  would 
be  [in  Daubechies’  notation] 


H-kr)  a  >  1,  j,k€Z, 


and  where  t  =  sample  time  interval. 


In  summary,  the  wavelets  are  basis  functions  for  the  wavelet  expansion  just  as  the  complex 
exponentials  (sinusoids)  are  the  basis  functions  for  the  Fourier  expansion.  The  continuous  wavelet 
expansion  of  a  function  (signal)  x(0  in  (91)  is  the  inner  product  defined  as 

/OO  .  _ 

x{t)n,C—^)dt 

■OO  a 


or  discretized^ 

W  [r(<)]  =  K.  =  !“!" *  r  -  kT)di  =  <  *,  > 

where  a  ->  oj,  T  ->  /:cir ,  and  >  1  (generally  2),  j,  k  e  2.  The  coefficient  outside  the  integral  is 
the  normalization  factor  such  that  the  L^-norm  i|V|l2  =1.  The  o^-term  is  the  magnification  where  _/ 
negative  and  large  corresponds  to  large  magnification,  and  consequently  aj  T  =  T/ai^  corresponds 
to  small  steps,  and,  therefore,  fine  details  are  discernible  (i.e.,  captures  high-fiequency,  short  tran¬ 
sient  characteristics). 


The  major  mathematical  difference  between  wavelets  and  the  Fourier  basis  functions  is  the 
requirement  that 

I  Tp{t)dt  =  0  , 

J  —  OO 


2.  Hoe  discretized  means  that  only  integral  shifts  in  scale  and  window  position  are  to  be  considered.  How¬ 
ever,  this  is  still  the  continuous  wavelet  transform  (CWT)  since  the  integral  and  Jt(f)  are  continuous  in  time, 
and  this  discretization  is  merely  for  analysis  ecotxxny. 
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which  allows  the  introduction  of  the  dilation  (or  scale)  parameter  in  order  to  make  the  time-frequen¬ 
cy  window  flexible  (Ref.  2).  However,  this  also  means  there  is  no  DC-component  in  wavelet  anal¬ 
yses,  but  then  again,  the  time-frequency  window  would  have  to  be  infinite  to  hold  the  DC- 
component.  What  this  finally  implies  is  that  the  final  expansion  can  only  reduce  the  signal  to  an  m- 
point  approximation  (for  an  m'*-order  wavelet)  in  time,  not  a  single  DC  component. 

The  wavelet  is  developed  not  from  just  a  single  mother  analyzing  function  as  in  the  Fourier 
case,  but  from  a  single  scaling  function  which  is  orthogonal  to  the  analyzing  wavelet  basis  function. 
In  contrast  to  the  STFT,  which  essentially  windowed  the  basis  functions  (or  transforms  only)  of  the 
expansion,  the  wavelet  truly  does  window  the  time  domain  and  the  Fourier  transform  directly.  This 
leads  naturally  to  a  filtering/decimation  process.  Consider  <)i(0  to  be  a  smoothing  function  which 
removes  irregularities  by  averaging  the  signal  locally  over  some  time  (or  sample)  interval.  If  one 
performs  the  inner  product  of  x{f)  and  \|r(0.  then  one  obtains  an  approximation  of  jt(0  within  the 
time  resolution  of  t(r(/),  i.e.,  <  xyf  >,  is  a  smoothed  version  of  x.  The  scaling  function  must  also 
cover  the  real  line  and  average  over  varying  intervals.  Consequently,  the  scaling  function  has  the 
same  form  as  the  wavelet,  i.e., 

4>j.k{t)=a~^(l>ia~H-kT)  a  >  1,  j,ke2, 

and  where  again  t  =  sample  time  interval.  However,  the  scaling  functions  only  window;  they  don't 
wiggle  or  oscillate. 

The  wavelet  function  extracts  the  detail  from  the  signal,  the  difference  between  the  original 
and  smoothed  versions,  such  that  the  signal  approximation  is  composed  of  the  smoothed  version 
plus  the  detail  removed  by  the  smoothing  operation.  Therefore,  two  function  spaces  are  established: 

<^GV'=»{0}c--  c1/2CKiCKoC  V_iCV_sC - ►Z,2(»)  =>  SCALE  SPACE 


and 


e  =>  {0}  ®  ®  1^2  ©  IVi  0  M^O  0  V/_1  ©  © - »  =►  DETAIL  SPACE 


These  spaces  are  related  in  that 


Vj^^  =  VjeWj, 
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i.e.,  the  next  finer  level  of  approximation  is  generated  by  the  orthogonal  sum  of  the  smoothed  ap¬ 
proximation  plus  its  detail  lost  due  to  smoothing  at  that  level.  Note  that  the  wavelet  spaces  Wj 
inherit  the  scaling  property,  i.e.,  jc(0  e  <=>  xiiJt)  s  from  Vy;  therefore,  x(t)  e  <=>  € 

Also  note  that  is  the  orthogonal  complement  of  Vj  in  V}.],  i.e.,  Wj ,  the  detail,  does  not  exist 
in  V^-;  it  exists  in  the  next  finer  level  of  approximation,  Consequently,  for  j  <  J, 

j-j-i 

Vj  =  V>  ®  0  Wj.t  ,  =  0lVi  , 

*=0  k£Z 


defines  a  multi-resolution  pyramidal  decomposition  where  all  of  the  spaces  are  orthogonal,  i.e., 

VjlVj',  WjlWy,  and  V^LW^. 

Note  that  Vj  is  decomposed  into  a  coarser  approximation  plus  its  lost  detail,  but  never  do  we 
decompose  the  detail  any  further,  as  shown  below.  This  way  the  energy  of  the  signal  is  preserved. 


WAVELET  PYRAMIDAL  DECOMPOSITION 
Since  the  interest  in  wavelets  is  for  the  on-line  analysis  of  a  signal  which  has  been  digitally 
sampled,  it  makes  sense  to  look  at  wavelets  from  a  filtering  (digital  signal  processing)  point  of  view. 
Let  a  =  2,  then  the  wavelet  will  have  binary  dilations  (V  nAt)  and  dyadic  translations  IkAt/Tf)  where 

Notice  that  dilating  a  signal  by  a  factor  of  2  is  equivalent  to  subsampling  in  the  sense  that  only 
half  the  samples  are  necessary  for  similar  resolution  in  a  given  time  interval  per  Nyquist.  Also 
notice  that,  because  of  the  duality  between  frequency  and  time,  for  a  binary  dilated  signal  (1/2  the 
frequency)  twice  the  time  interval  is  required  to  sample  that  signal  as  compared  to  the  original 
signal,  non-dilated  version,  for  a  similar  level  of  resolution.  Consequently,  one  can  lowpass  filter, 
determine  the  detail  lost  by  subtracting  the  lowpassed  version  from  the  original,  and  then  decimate 
by  a  factor  of  2  and  accumulate  more  samples  (look  at  longer  time  intervals)  to  get  any  resolution 
in  frequency  as  desired  at  the  expense  of  time  resolution.  This  is  the  MRA  described  above. 
However,  this  MRA  requires  the  two  basis  functions  described  by  the  wavelet  approach.  In  this 
case,  ^  performs  the  lowpass  filtering  and  decimation  to  coarsen  the  approximation  (smooth  the 
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signal),  and  y  simultaneously  performs  the  highpass  filtering  at  the  previous  level  to  preserve  the 
detail  (fluctuations)  which  will  be  lost  going  to  this  coarser  level  of  approximation.  The  wavelet 
analysis  then  processes  the  data  according  to  the  grid  as  shown  below. 

f 

k 

Original  Samples 


I 


Analysis 


Synthesis 


t 


WAVELET  CALCULATION  GRID 

Note  again  that  no  E>C  processing  is  performed,  that  the  grid  is  much  more  sparse  than  that  for  the 
STFT,  and  finally  that  not  all  frequencies  are  processed:  rather,  the  frequency  axis  is  scaled  accord¬ 
ing  to  a  base  2  logarithm  due  to  the  half-band  filtering  and  subsequent  decimation  by  two. 

As  mentioned  above,  the  wavelet  (detail)  function  spaces  and  the  scale  function  spaces  are 
related.  As  filters,  this  places  some  unique  requirements/restrictions  on  the  coefficients  such  that 
perfect  reconstruction  is  possible.  These  restrictions  result  in  a  class  of  filters  known  as  quadrature 
mirror  filters  (QMF).  Quadrature  mirror  filters  (see  schematic)  are  filters  G  and  H  which,  for  all 
signals  X  of  finite  energy,  produce 


liVL|p  +  l|y„||'=||;f  IP 

where  is  the  output  from  the  lowpass  filter  C  and  is  the  output  from  the  highpass  filter  ff, 
bo  A  of  which  are  operators  which  map  l\z)  =»  /^(22),  where  /  is  the  Hilbert  space  such  that 
^  |jc(n)  I  <  oofor  nE2.  In  the  case  of  perfect  halfband  filters,  it  is  obvious  that  the  coarser 
Spi^ximation  contains  only  those  frequencies  below  nfl  of  the  original  signal  jc(n),  while  the 
highpassed  signal  yu  contains  only  those  fiequencies  above  iU2.  Consequently,  decimation  is 
justified  for  both  the  coarser  approximation  and  the  detail,  holding  the  number  of  coefficients 
({A] },  {Dg})  constant  at  the  number  of  original  data  samples.  Note  that  reconstruction  is  perfect, 
and  this  may  be  expressed  by 
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I  =  G'G  +  H'H 

where  G'  and  if  are  the  adjoints  of  the  filters  G  and  H,  respectively,  mapping  l\2Z)  =>  In 
addition,  note  that  if  g(/i)  denotes  the  impulse  response  of  the  lowpass  filter  with  frequency 
response  Gfco),  a  highpass  filter  can  be  obtained  by  translating  G((i})  by  it-radians  (i.e.,  replacing  o> 
by  ©  -  Ti)  and  correcting  for  phase.  Therefore,  the  highpass  frequency  response  //(©)  =  G*(ffl  -  Tt), 
or  Mji)  =  e'^  g(l  -  n)  =  (-l)"g(l  -  n),  which  is  true  for  finite  impulse  response  (FIR)  filters.  The 
following  block  diagram  illustrates  the  quadrature  minor  filtering  process. 


QUADRATURE  MIRROR  FIITERING 

Daubechies  (Ref.  5)  has  shown  that  such  filters  can  be  constructed  with  all  of  the  desirable 
characteristics,  i.e.,  compact  support  in  both  the  frequency  and  time  domains,  orthogonality,  mini¬ 
mal  coefficients  (no  more  than  the  number  of  original  samples)  and  perfect  reconstruction.  In  this 
case  the  scale  functions  and  wavelet  functions,  respectively,  are 

(^ij.*(7i)  =  -  k) 

k 

and 

k 

which  are  known  as  two  scale  difference  equations  or  dilation  equations,  and  h(n)  =  (-1)"  g(n  - 1). 

The  construction  of  wavelets  then  begins  with  the  scaling  function  Newland  (Ref.  8,  pp. 
308-321)  provides  an  excellent  explanation  as  to  how  these  filter  coefficients  are  determined  based 
upon; 
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•  CONSERVATION  of  AREA  (Unique  Scaling  fVinction): 

iV-l 

=  V2 

t:=0 

•  ACCURACY: 

N-l  „ 

=  m  =  0:1,2,...,:^ -1 

i=Q  ^ 

•  and  ORTHOGONALITY: 

iV-l  N-l 

X]  9i9k+2m  =  0,  and  =  1, 

*=o  i=0 

where  N  is  the  number  of  coefficients  desired.  These  conditions  result  in  a  set  of  N  +  1  simultaneous 
equations  in  terms  of  and  finally,  =  (- 1  )*  gt-  Assuming  that  such  filters  are  available,  the  trans¬ 
formation  of  time  domain  data  to  a  time-scale  representation  resulting  from  convolutions  of  the  sig¬ 
nal  with  these  filters  follows  the  pyramidal  scheme  as  coded  by  Press  (Ref,  9)  and  illustrated  below. 


DISCRETE  WAVELET  TRANSFORM  PYRAMIDAL  ALGORITHM 


The  most  unique  aspect  of  wavelet  analysis  is  that  the  time-domain  data  is  transformed  not 
into  the  familiar  2D  amplitude-frequency  plane  as  with  Fourier  analysis,  but  rather  into  a  3D 
amplitude-scale-time  domain.  Again,  scale  refers  to  the  broadening  of  the  basis  function  to  fill  the 
time  (or  space)  interval  associated  with  the  filter  order,  i.e.,  an  m^-order  wavelet  filter  requires  m 
data  points  for  the  convolution  process  regardless  of  the  level  of  the  analysis.  Consequently,  the 
basis  function  is  more  compact  for  the  first  level  of  wavelet  analysis  and  broadens  in  time  (or  space) 
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with  each  successive  level  of  the  analysis  due  to  the  decimation  process  which  lengthens  the  time 
(or  space)  interval  between  successive  data  points.  The  inverse  of  this  scale,  or  period-of-existence. 
of  the  basis  function  takes  on  a  similarity  \a  frequency  in  a  mechanical  sense  if  the  basis  function 
has  a  single  zero-crossing  as  in  the  Haar  and  4-tap  Daubechies  wavelets.  This  is  due  to  the  fact  that 
a  material  fiber  would  undergo  a  tensile/compressive  cycle  at  a  rate  of  once  per  scale  of  the  basis 
function,  similar  to  that  of  a  sine  wave  in  Fourier  analysis. 

This  can  be  better  shown  if  one  considers  a  sine  wave  signal  composed  of  two  harmonics  plus 
the  fundamental  for  Fourier  analysis  sampled  at  8  times  per  cycle  for  the  highest  harmonic,  and  the 
same  except  a  square  wave  for  Haar  wavelet  analysis  as  shown  in  Fig.  6.  The  Fourier  analysis  ap¬ 
plied  to  the  sine  waves  illustrates  each  of  the  harmonics  directly  in  the  Fourier  (amplitude  vs. 
frequency)  plane  as  coefficients  of  the  sine  waves  making  up  the  signal  represented  as  spikes  of  the 
appropriate  amplitude  at  the  appropriate  frequency.  The  wavelet  analysis  applied  to  the  square 
waves  illustrates  these  harmonics  in  the  3D  time-scale  plane  as  detail  coefficients,  which  are 
constant  valued  functions  of  the  appropriate  amplitude  at  the  appropriate  scale  of  the  basis  function 
extending  for  all  time  (or  space)  of  the  sample  block.  These  functions  are  constant-valued  in  time 
since  the  square  wave  signal  is  periodic,  just  as  the  Fourier  coefficients  are  constant  in  time  (or 
space)  since  Fourier  analysis  assumes  stationary  signals  and  all  temporal  (or  spatial)  dependence  is 
lost  following  the  FFT. 

Recall  that  the  wavelet  detail  coefficients  are  the  amplitudes  of  the  fluctuating  continuous 
wavelet  basis  functions,  just  as  the  FFT  coefficients  are  the  amplitudes  of  the  continuous  sine  wave 
basis  functions,  to  be  combined  to  synthesize  the  respective  signals.  Similarly,  the  wavelet  scale 
coefficients  are  the  decimated  signal  samples,  i.e.,  a  smoothed  representation  of  the  signal  at  a 
coarser  scale.  For  the  Haar  wavelet  applied  to  a  square  wave  oversampled  at  8  samples  per  cycle, 
the  wavelet  amplitudes  are  all  zero  except  the  third  level  of  decomposition  (for  512  samples,  this  is 
level  6  of  the  9  possible,  i.e.,  2^  =  5 12,  level  9  being  the  original  samples),  which  is  constant  for  all 
values  of  time  (or  space)  for  the  periodic  square  wave.  Consequently,  the  Haar  wavelet  does  for 
square  waves  what  Fourier  analysis  does  for  sine  waves,  i.e.,  determines  the  exact  scale  or 
frequency  match,  respectively,  for  a  signal  composed  of  the  respective  basis  functions.  The  Haar 
wavelets,  of  course,  can  only  provide  such  perfect  scale  matches  if  the  data  is  properly  aligned  to 
harmonics  and  sampled  appropriately,  but  this,  too,  is  no  different  than  the  FFT,  where  pure  spikes 
will  occur  only  if  the  data  has  integral  frequencies  and  is  sampled  with  zero  phase.  Otherwise, 
smearing  occurs  in  the  FFT  spectrum,  and  the  wavelet  analog  is  to  have  numerous  non-zero,  and 
typically  oscillatory,  details  at  scales  below  the  highest  harmonic  in  the  signal.  It  should  be  noted 
that  at  the  third  level  of  decomposition,  the  wavelet  analysis  has  decimated  the  original  signal  by  a 
factor  of  8,  the  exact  factor  of  the  samples  per  cycle,  leaving  only  one  sample  per  cycle.  This  holds 
for  other  sampling  factors  as  well.  In  other  words,  the  log2  of  the  number  of  samples  per  cycle  is  the 
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decomposition  level  which  will  be  non-zero  if  the  signal  is  a  series^  of  the  basis  function  of  the 
wavelet.  At  this  level,  the  reciprocal  of  the  interval  between  the  remaining  samples  is  the  frequency 
of  the  square  wave. 

It  should  also  be  pointed  out  that  the  choice  of  the  wavelet  basis  function  is  equivalent  to 
making  an  assumption  about  the  composition  of  the  signal.  The  use  of  a  Haar  wavelet  for  a  square 
wave  is  making  the  obvious  assumption  that  a  square  wave  is  made  up  of  the  Haar  basis,  just  as  in 
the  case  of  using  the  FFT  makes  the  assumption  that  the  signal  is  composed  of  sine  functions.  Each 
basis  applied  to  the  wrong  signal  will  generate  broad  spectra  which  tend  to  yield  little  or  no 
information  about  the  signal  composition.  The  potential  of  wavelets  is  that  we  now  have  a  broad 
selection  of  basis  functions  to  apply  to  a  problem  to  detect  characteristics  which  are  smeared 
utilizing  only  the  FFT,  or  to  perform  pattern  recognition  analysis  where  we  generate  a  new  basis 
dependent  upon  the  pattern  of  interest. 

43  APPLICATION  TO  MODEL  SIGNAL 

Consider  the  same  signal  used  for  the  STFT  analysis  discussion  of  Section  3.0  (see  Fig.  1). 
Application  of  the  Daubechies  four-term  wavelet  (Daub4)  to  this  signal  results  in  Fig,  7.  Again,  the 
original  signal  (wavelet  level  10  since  the  number  of  samples  =  1024  =  2*^)  is  shown  at  scale  index 
2  and  the  first  wavelet  (level  9)  is  shown  at  scale  index  3.  Therefore,  subtracting  2  from  the  scale 
index  in  these  3D  wavelet  plots  will  yield  the  wavelet  iteration  number,  which,  when  subtracted 
from  10  (the  highest  wavelet  level  for  these  data)  yields  the  wavelet  scale.  Now  the  signal  has  three 
distinct  frequencies,  /q  =10  Hz,  /]  =  80  Hz,  and  /2  =  128  Hz.  The  sampling  rate  was  2.56*F„^  or 
327.68  Hz.  The  Itrst  wavelet  filters  the  signal  and  subsamples  by  a  factor  of  2,  which  means  at  level 
9  the  sample  rate  of  the  wavelet  scale  is  effectively  163.84  Hz.  This  is  essentially  twice  /],  and 
therefore  the  80-Hz  signal  is  shown  as  a  sinusoid  at  scale  index  4,  which  is  level  8  (this  was  veriiied 
using  a  single  80-Hz  component  sampled  at  128/80  ■  2.56  or  4.096  samples  per  cycle).  The  trace  at 
scale  index  3  (level  9)  should  be  the  /^-modulated  128-Hz  cosine  wave.  Recall,  the  wavelet  shows 
the  detail  lost  by  the  filtering  at  the  previous  level.  The  interpretation  of  this  trace  is  not  fully  un¬ 
derstood  by  the  authors  and  it  is  suspected  that  this  trace  is  noisy  due  to  the  wavelet  basis/signal 
mismatch.  The  Gaussian  impulses  are  observed  as  spikes  at  sample  time  indices  320  and  576  which 
broaden  with  successive  wavelet  iterations.  Note  that  the  temporal  characteristics  were  preserved 
with  the  wavelet  transform,  unlike  the  FFT,  and  the  spikes  occur  at  the  peaks  of  the  impulse  even 
though  the  impulses  were  rather  broad  (64  samples  each).  This  is  due  to  the  siope  change  from  one 
side  of  the  peak  to  the  other,  and  therefore  an  edge  is  detected  by  the  wavelet  at  this  point.  The  trace 
at  scale  index  6  (level  6)  is  an  artifact  of  sampling  the  10-Hz  signal  at  a  factor  of  32.768  (2.56) 


3 .  Here  series  means  a  series  in  time  (or  space)  of  wavelet  basis  functions  laid  end  to  end. 
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single  10- Hz  component  sampled  at  32.768  samples  per  cycle).  This  artifact  appears  to  be  due  to  a 
beating  between  the  wavelet  zeros  and  the  signal  zeros  at  this  level,  but  this  could  not  be  verified 
by  the  author.  This  artifact  persists  regardless  of  the  order  of  the  wavelet,  from  Haar  to  Daub20.  The 
only  significant  changes  going  to  other  wavelet  orders  were  the  coarsening  of  the  traces  with  the 
Haar  wavelet  and  the  smoothing  of  the  traces  with  the  Daub20  wavelet,  and  the  more  pronounced 
spikes  at  the  Gaussian  impulse  positions  due  to  less  noise  with  the  Daub20.  The  Daub20  wavelet 
results  are  shown  in  Fig.  8.  Note  that  neither  Fig.  7  nor  Fig.  8  contains  any  information  about  the 
amplitude  modulation  aspect  of  the  signal.  This  is  a  deficiency  of  the  wavelet  analysis,  at  least  as  it 
is  presently  understood  and  interpreted  by  the  authors. 

The  STFT  provided  results  which  were  considered  more  familiar  to  that  of  the  wavelet 
analysis  in  interpretation  aspects,  but  the  wavelet  analysis  is  much  faster,  has  less  impact  on  the 
hardware  requirements  for  sampling,  and  provides  suitable  information  once  the  interpretation  is 
worked  out.  In  addition,  the  wavelet  could  be  used  in  multi-resolution  schemes  to  properly  filter  the 
signal  for  FFT  analysis  or  to  complement  the  FFT  by  searching  for  impulses  or  other  discrete  non¬ 
stationary  events.  Further  research  into  the  interpretation  of  wavelet  data  must  be  provided,  and 
perhaps  the  new  harmonic  wavelet  of  Newland  (Ref.  8)  will  shed  some  light  in  this  area, 

4A  APPLICATION  TO  PULSE-ECHO  ANALYSIS 


Signal  discontinuities  are  prevalent  in  real  signals,  and  one  such  signal  which  relies  on  discon¬ 
tinuities  is  that  of  Pulse-Echo,  which  is  a  technique  for  measuring  distances  in  materials  nonintru- 
sively  through  the  use  of  sound  waves.  In  brief,  a  sound  impulse  is  propagated  into  the  material 
through  an  appropriate  acoustic  coupling  and  a  sensor  detects  the  echo  or  reflected  wave  caused  by 
the  abrupt  change  in  the  index  of  refraction  at  the  opposite  edge  of  the  material  (this  can  be  the 
interface  between  two  different  materials  or  an  air/material  interface).  The  thickness  d  of  the  mate¬ 
rial  is  then  found  as  * 

d  =  c- 


where  c  is  the  characteristic  speed  of  sound  for  the  material,  and  t  is  the  round  trip  time  from 
initiation  of  the  pulse  through  the  material  to  the  discontinuity  in  the  index  of  refraction  and  back 
through  the  nuterial  to  the  sensor.  The  sensors  are  typically  piezoelectric  crystals  of  appropriate 
natural  frequency  which  act  as  the  transmitter  when  electrically  pulsed  and  act  as  the  receiver  when 
stimulated  by  the  returning  sound  wave  echo  (hence,  the  name  Pulse-Echo).  The  signals  resemble 
an  impulse  with  ringing,  (similar  to  a  damaged  bearing  signal  (Ref.  4)  and  were  simulated  by  single 
sided  sine  functions  with  a  reduced  amplitude  for  the  echo  to  simulate  material  transmittance 
properties.  The  functions  used  were  of  the  form 
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a;[t]  =  u_i(i  -  ti)j4i 


sin{i  —  Ti) 

(<  -n) 


+  u_i(i  -  ti)A2 


sin((  -  Ti) 


where  t]  was  near  the  beginning  of  the  ensemble,  ^  midpoint  of  the  ensemble,  and  Ai 

was  0.5A| .  The  first  sine  function  was  spread  over  the  whole  interval  once  started,  while  the  second 
sine  function  started  such  as  to  peak  at  the  exact  midpoint  and  build  upon  the  first  impulse,  hence 
the  unit  step  functions  -  T).  The  signal  was  sampled  at  16  samples  per  cycle  for  the  1-Hz  sine 
wave  component  of  the  sine  functions,  and  1024  point  ensembles  were  used.  Random  noise  was 
added  to  the  signal,  and  various  noise  floor  amplitudes  were  evaluated.  The  closeness  in  time  be¬ 
tween  the  simulated  pulses  could  be  controlled  to  simulate  varying  material  thickness.  The  test  pa¬ 
rameters  are  summarized  in  Table  1,  and  the  wavelet  3D  detail  plots  are  shown  in  Fig.  9. 


In  pulse-echo  measurements,  assuming  the  material  speed  of  sound  characteristics  are  well 
known,  the  most  important  parameter  is  an  accurate  measurement  between  the  main  pulse  and  it’s 
echo.  Any  point  on  the  pulses  can  be  used  as  a  reference  point,  since  the  same  sensor  is  used  for 
transmitting  and  receiving  it  will  have  the  same  response  function,  but  it  is  necessary  to  be 
consistent.  Since  real  signals  generally  have  a  significant  noise  floor,  it  is  usually  most  accurate  to 
use  the  peaks  in  the  pulses  as  the  reference  points.  However,  most  pulse-echo  instruments  use  a 
threshold  test  and  capitalize  on  the  rapid  decay  of  the  response  because  detecting  the  true  peak 
location  is  difficult  electronically.  The  use  of  wavelets  may  change  this  situation  due  to  their  ability 
to  detect  edges  and  rapid  slope  changes  in  a  signal. 


Thble  1.  Pulse-Echo  Test  Cases 
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In  pulse-echo  measurements,  assuming  the  material  speed  of  sound  characteristics  are  well 
known,  the  most  important  parameter  is  an  accurate  measurement  of  the  time  interval  between  the 
main  pulse  and  its  echo.  Any  point  on  the  pulses  can  be  used  as  a  reference  point,  since  the  same 
sensor  is  used  for  transmitting  and  receiving  it  will  have  the  same  response  function,  but  it  is 
necessary  to  be  consistent.  Since  real  signals  generally  have  a  significant  noise  floor,  it  is  usually 
most  accurate  to  use  the  peaks  in  the  pulses  as  the  reference  points.  However,  most  pulse-echo 
instruments  use  a  threshold  test  and  capitalize  on  the  rapid  decay  of  the  response  because  detecting 
the  true  peak  location  is  difficult  electronically.  The  use  of  wavelets  may  change  this  situation  due 
to  their  ability  to  detect  edges  and  rapid  slope  changes  in  a  signal. 
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Various  Daubechies  wavelets  were  applied  to  the  signal  with  the  intent  of  locating  the  peaks  as 
consistent  reference  points  for  determining  the  time  interval.  The  Haar  (DaubZ)  wavelet  provided  the 
best  results,  and  was  used  for  this  analysis.  The  FFT  applied  to  the  signal  of  Fig.  9a  merely  shows  a 
pulse  near  DC  and  offers  no  timing  information  (see  Pig.  10).  It  does  illustrate  that  the  signal  is 
appropriately  sampled,  howevar,  since  the  spectrum  is  essentially  zero  well  before  the  Nyquist 
foldover  frequency.  Figure  1 1  illustrates  the  ability  of  the  wavelet  to  detect  the  location  of  the  peaks 
at  wavelet  levels  4  and  3.  Level  4  is  significant  since  at  this  level  the  signal  is  sampled  at  essentially 
the  Nyquist  criteria.  Even  in  the  presence  of  significant  noise  (Fig.  9c),  the  Haar  wavelet  was  able  to 
detect  the  peaks  in  either  of  these  two  levels.  While  the  fourth  level  has  a  second  peak  which  becomes 
more  pronounced  with  increasing  noise,  the  third  level  clearly  distinguishes  which  peak  in  level  4  is 
the  appropriate  peak  of  the  pulse  for  timing  purposes.  However,  the  ability  of  wavelets  to  determine 
the  peak  locations  is  tied  to  whether  or  not  the  peak  was  sampled  and  if  the  sample  index  n  is  a  factor 
of  2f,  This  is  illustrated  in  Fig.  12a,  where  the  starting  point  was  at  n  =  169,  and  a  maximum  error  of 
13  percent  (for  level  3;  level  4  is  7  percent)  is  incurred  due  to  the  subsampling  resulting  in  the  peak 
being  missed  slightly.  This  is  equivalent  to  the  scalloping  effect  of  the  FFT  when  the  signal  frequency 
does  not  correspond  to  the  FFT  bin  center  frequency  but  rather  is  offset  from  the  center  of  the  bin. 
Figure  12b  takes  this  situation  to  the  extreme  by  narrowing  the  distances  between  the  pulses  and 
having  n  a  non-integral  power  of  two.  This  is  probably  the  minimum  spacing  that  would  be  consid¬ 
ered  a  typical  pulse-echo  signal  (for  thinner  materials,  higher  frequencies  (tighter  pulses)  are 
typically  used  to  space  out  the  pulses  for  better  time  resolution  of  the  peaks).  The  peaks  are  still  well 
defined  in  level  4,  while  level  3  has  rounded  and  smeared  the  peaks. 

In  summary,  the  wavelet  analysis  provided  a  consistent  set  of  timing  points  for  improving  the 
resolution  of  the  pulse-echo  measurement  technique,  especially  when  the  signal  is  buried  in 
significant  noise.  The  limitation  of  the  approach  appears  to  be  the  pulse  peak  location  relative  to  an 
integral  power  of  two  sample  points  within  the  ensemble.  This  could  be  resolved  by  performing  a 
sliding  wavelet  analysis  such  that  the  peak  would  become  more  pronounced  as  it  fell  on  an  integral 
power  of  two  sample,  and  then  timing  between  two  different  wavelets  at  the  same  level.  Work  is 
being  performed  to  make  the  sliding  wavelet  computationally  efficient,  and  to  calculate  a  particular 
level  without  having  to  calculate  all  previous  levels  serially. 

5,0  MALVAR  WAVELET  TRANSFORM 

Traditional  signal  analysis  using  the  Fourier  transform  relies  heavily  on  the  fact  that  the  data 
under  study  is  stationary  in  some  sense.  The  usual  display  of  the  spectra  uses  a  single  window  and 
provides  no  way  to  show  variations  in  the  spectra  within  the  window;  usually  only  one  window  size 
is  used  throughout  the  analysis,  tacitly  assuming  staiionarity  at  least  on  the  order  of  several  window 
sizes.  Window  effects  such  as  leakage  are  controlled  by  reducing  window  influence  at  the  ends  of 
the  window  interval,  thus  also  discarding  much  of  the  influence  of  the  signal  data  at  these  transition 
points  and  making  reconstruction  problematic.  Noise  reduction  techniques  are  based  on  averaging 
over  several  windows  which  are  considered  to  be  stationary  in  a  statistical  sense.  Signals 
representing  non-stationary  phenomena  must  be  treated  so  as  not  to  discard  any  parts  since  it  cannot 
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be  assumed  that  “one  part  is  much  like  another.”  This  will  require  basis  functions  which  vary  in 
time  locality  as  well  as  frequency  range.  These  basis  functions  both  within  and  between  windows 
must  allow  complete  and  stable  reconstruction  from  the  analysis  coefficients. 

Malvar  wavelets  are  a  natural  development  from  windowed  Fourier  transforms.  The  signal 
sit)  is  first  sliced  up  in  time  using  window  functions  Wj<r)  to  obtain  a  series  of  windowed  signals 
Each  of  the  windowed  functions  is  then  subjected  to  Fourier  analysis.  This  closely  resem¬ 
bles  the  Gabor  (Ref.  10)  approach,  since  the  Fourier  transform  leads  to  calculating  fs  <  s.  Wj^ 
>with  the  wavelet  Wfjit)  =  The  windows  wjit)  are  chosen  to  give  several  advantages. 

Gabor  used  a  single  Gaussian  with  offsets  for  Wj{t)  which  has  a  number  of  disadvantages.  The 
Gaussian  does  not  have  compact  support  and  the  algorithms  must  be  designed  with  error  analysis 
in  mind.  But  more  serious  is  the  problem  of  accounting  for  the  overlap  of  the  windows;  this 
becomes  clear  if  we  try  to  reconstruct  sfr)  from  the  spectra.  The  phenomena  known  as  “leakage” 
in  Fourier  analysis  are  caused  by  windowing  which  does  not  completely  include  a  feature  (such  as 
a  complete  cycle)  in  a  single  window,  but  instead  divides  the  influence  across  two  adjacent  win¬ 
dows.  What  is  needed  for  reconstruction  (and,  as  we  will  see,  for  control  of  some  other  influences) 
is  that  we  have  orthogonality  between  the  windows  as  well  as  within  them. 

Malvar  studied  lapped  Fourier  transforms  (Ref.  1 1),  which  can  be  regarded  as  overlapping 
rectangular  windows,  with  the  intent  of  bringing  about  orthogonality  between  windows.  Actually, 
rectangular  non-overlapping  windows  are  already  orthogonal.  Newland  exploited  this  fact  to  use 
the  Fourier  transform  of  a  rectangular  window  to  construct  a  family  of  wavelets  based  on  the  sync 
function  (Ref.  S).  But  the  edge  effects  of  the  rectangular  window  are  severe.  Figure  13  shows  a 
reconstruction  of  an  acoustic  signal  from  spectra  given  in  several  adjacent  windows  after  discarding 
several  of  the  least  significant  spectral  coefficients;  note  the  cusp  at  several  interval  boundaries 
caused  by  the  slow  roll-off  of  the  Fourier  spectrum  of  the  rectangular  window.  Malvar  found 
orthogonality  conditions  on  the  lapped  transforms  which  can  be  described  similarly  to  a  windowed 
Fourier  transform.  But  his  window  extends  outside  the  “primary”  interval  and  overlaps  the  adjacent 
window.  Figure  14  shows  how  the  windows  of  two  adjacent  intervals  are  related  while  Fig.  15 
shows  what  a  typical  Fourier  basis  function  looks  like  after  being  windowed.  The  smooth  edges  on 
the  windows  are  reminiscent  of  window  “carpentry”  which  has  been  used  to  reduce  leakage  and 
other  window  effects  (Ref.  12);  but  the  older  techniques  did  not  overlap  the  windows  nor  account 
for  inter-window  signal  effects. 

Coifman  and  Meyer  (Ref.  13)  modified  Malvafs  approach  to  obtain  variable  width  windows 
while  continuing  to  maintain  orthogonality  within  and  between  the  windows.  Their  approach  still 
retains  the  faster  roll-off  in  window  effects  obtained  by  using  smooth  window  edges. 

Only  the  window  shapes  of  part  of  a  possible  basis  are  shown  in  Fig.  16.  This  basis  would  be 
suitable  to  represent  a  signal  which  was  almost  stationary  in  region  1  and  4  while  significant 
transient  behavior  is  captured  in  regions  2  and  3.  We  must  emphasize  again  that  the  transform 
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transient  behavior  is  captured  in  regions  2  and  3.  We  must  emphasize  again  that  the  transform 
coefficients  on  each  of  the  intervals  do  not  interact;  the  orthogonality  ensures  this.  Furthermore,  die 
effect  of  a  feature  extending  between  two  adjacent  windows  is  captured  such  that  it  can  be 
accurately  reconstructed  in  spite  of  round-off  or  truncation  errors. 

We  illustrate  the  construction  of  the  Meyer-Coifinan-Malvar  discrete  transformation  using  a 
cosine  type  IV  basis  on  each  basic  window  interval.  Other  bases,  including  the  Fourier  exponential 
basis,  can  be  used,  but  the  symmetries  needed  to  ensure  the  behavior  referred  to  in  the  previous 
paragraph  are  only  simple  in  the  cosine  type  FV  case. 

Figure  17  shows  the  nomenclature  for  both  the  continuous  and  the  discrete  versions  of  the  Malvar 
wavelets.  The  Cj  represent  the  “classic”  interval  end-points  for  non-overlapped  windows.  The  jth 
interval  extends  from  cj  to  cj+i  and  has  length  Ij. 

The  window  w'W  corresponding  to  thejth  interval  is  constructed  by  modifying  the  rectangular 
window  using  the  bell  profile  function  p/(/)  shown  in  Fig.  18.  The  result  is  that  the  window  n/fr) 
extends  Oj  to  the  left  of  cj  and  ocy^ito  the  right  of  c^\.  The  window  is  also  modified  within  the 
original  basic  interval  as  shown  in  the  Fig.  18.  The  “mirror  image”  of  P^(0  is  used  to  modify  the 
right  end  of  the  (/  -  l)th  interval.  The  particular  of  Fig.  18  ensures  \J{i)  has  continuous 
derivatives  even  at  cj  ±  ay  P'(r)  can  be  chosen  to  make  as  smooth  as  desired. 

Discretization  is  simplified  by  rescaling  the  entire  problem,  if  necessary,  so  that  the  sample 
points  on  thejth  interval  occur  at  integers  +  n,  0  jS  n  <  Ij  and  the  Cj  occur  exactly  halfway 
between  Uj^jand  Oj.  Rj  is  the  largest  integer  less  than  aj  1/2. 

p^(/)  is  discretized  so  that  Pi(  =  +  n)  and  other  functions  of  /  are  similarly  characterized. 

The  discrete  window  function  wi  for  the  Jth  interval  which  results  is  given  by 


^  _  i  1  if  n  6  [flj.fi  -  1  -  flj+i], 

"  ifnG[/j-flj+i,/j4-flj+i-l]. 

'  0  otherwise. 


The  basis  functions  for  conventional  discrete  cosine  type  IV  transforms  are  given  in  continu¬ 
ous  and  discrete  form  by 

v4(0  =  J]'  ;  e  *  e  Af, 0  <  f  <  fj 


A,n  =  y  0  <  <  h 
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Thus,  the  continuous  and  discrete  basis  function  are 

=  <<Pk,n 

This  definition  does  not  yet  define  t]//  „  off  the  interval  from  Cj  to  cj+i,  but  the  symmetry  of  the 
basis  function  tpJlit)  at  the  end  points  and  a  symmetry  requirement  on  the  ^^(0  completes  the 
definition. 

=  V^.ni  'fii.Vj+n  =  ~‘Pk,n 

+  (j9in-i)"  =  1  ;  G  2,  n  6  [-/?> .  Rj  -  1] 

The  synunetry  in  Pj^  above  is  exactly  equivalent  to  the  conditions  obtained  by  Malvar  which 
ensure  orthogonality  between  basis  functions  in  two  adjacent  windows.  The  orthogonality  of  the 
discrete  basis  functions  within  and  between  windows  can  be  proved  exactly. 

The  discrete  basis  function  is  given  directly  in  terms  of  P;(  by 

ifn6  [--Rj./Zj  -  1], 

.  if  "  €  [ij  -  Rj,l}  +  R)  —  1]. 

Now  in  order  to  perform  the  transform  of  a  signal  s„,  we  need  the  transform  coefficients 
B )  ■  express  the  scalar  product  in  terms  of  a  new  signal  S  as  follows; 

V  ;  +  n  / 

■St  = 

nis— ilj 
Bj-l 

n=0 

ij-Wj+i-l 

Vj-l 

n=ij--Rj+j 

'.-1 

n=0 
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Thus,  using  the  new  signal  the  transform  becomes 


Sj[  —  {*o,+ni  IPfc.n) 


which  is  an  ordinary  discrete  cosine  transform  on  basic  interval  j.  The  special  signal  is  said  to 
be  obtained  by  “folding”  the  signal  at  cj  and  cy  ^  |.  As  expected,  the  transform  for  one  window 
includes  the  influence  of  part  of  the  signal  in  the  adjacent  windows;  it  obtains  a  weighted  average 
behavior  at  the  end-points  which  smooths  the  transition  across  the  window  edges.  A  similar 
“unfolding”  rule  is  obtained  easily  and  blends  together  the  transforms  on  two  adjacent  windows  to 
reconstruct  the  signal  without  artifacts. 

In  summary:  To  perform  the  Malvar  transform,  fold  at  each  af. 

(  +  Soj-n-l^n-1  if  ”  G  [0,  Rj  —  1], 

ioj+n  =  <  if  n  G  [-/ey,-!]. 

I  Sa^+n  otherwise. 

then  Discrete  Cosine  Transform  ^  on  each  interval  [oj,  oj+ij-  !]• 

To  invert  the  Malvar  transform,  first  DCT  as  before  (the  DCT-IV  being  used  here  is  its  own 
inverse),  then  unfold  at  each  Oj. 


Sflj+n 


^aj+n^  —  Soj-n-tM-n-l  if  "  €  [O,  Rj  —  1], 
"i"  *aj+n^n-l  if  n  G  \—Rj,  —1]. 
So,+n  otherwise. 


5.1  THE  WAVELET  UBRARY 

Since  we  may  use  windows  of  rather  arbitrary  size  in  the  signal,  it  is  clear  that  several  choices 
of  the  basis  intervals  and,  consequently,  of  the  basis  itself  are  possible.  The  collection  of  all  the 
Malvar  wavelets  on  all  the  intervals  in  one  arrangement  of  basis  intervals  is  a  basis  for  the  I? 
functions  on  the  real  line.  But  there  are  many  choices  for  the  basis  intervals,  so  there  are  many  basis 
sets  for  The  collection  of  all  these  “bases”  is  called  a  wavelet  library.  The  question  is:  Which 
basis  from  the  library  is  the  best  one? 

To  see  how  to  compare  two  bases  from  the  library,  we  consider  the  simple  case  represented  in 
Fig.  19,  where  we  have  drawn  only  the  windows  for  one  basis  rising  above  the  axis  and  a  window 
for  the  second  basis  but  reflected  below  the  axis.  If  two  bases  differ  only  in  having  the  two  upper 
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windows  in  one  replaced  by  the  lower  window  in  the  other,  then  we  need  only  decide  whether  it  is 
bener  to  use  the  single  window  (thereby  assuring  the  signal  is  nearly  stationary)  or  to  use  the  two 
smaller  windows  (thereby  more  accurately  capturing  any  transient  behavior). 

Wesfned  and  Wickerhauser  (Ref.  14)  have  applied  this  idea  in  acoustic  signal  processing  to 
recognize  speech  signals.  They  measured  information  in  transform  for  the  two  intervals  together 
against  that  in  the  single  interval  transform  with  the  formula  for  “spectral  entropy”  given  by 

//(Si)  =  -Y^Pilogpk  where  p*  =  (5i)^/|j:^||^. 

This  is  one  of  several  information  measures  which  can  be  used  to  “split  or  merge”  the 
windows  to  obtain  a  minimal  entropy  (maximal  information)  basis  for  the  signal. 

An  effective  and  fast  implementation  of  the  idea  is  to  first  obtain  Malvar  transfoims  on  rather 
small  windows,  say  Ij  =  32.  Then  measure  the  entropy  on  each  pair  of  intervals  separately  and 
grouped  as  one  with  Ij  =  64.  If  the  entropy  decreases,  keep  the  larger  and  discard  the  smaller. 
Continue  grouping  the  larger  intervals  in  pairs  and  replacing  if  the  entropy  decreases.  The  entropy 
measure  and  others  which  are  suitable  will  cause  this  process  to  obtain  the  optimal  basis  in  time 
comparable  to  that  required  by  the  discrete  cosine  transform.  Also  the  transform  coefficients  for  the 
two  smaller  windows  are  used  recursively  to  efficiently  compute  the  transform  in  the  larger 
window. 

Figure  20  shows  a  speech  signal  reconstructed  from  an  optimized  basis.  The  vertical  lines  rep¬ 
resent  the  optimal  windows  for  minimum  entropy.  Note  that  the  more  nearly  periodic  the  signal  is, 
the  longer  the  window.  This  is  exactly  the  behavior  we  hope  to  establish  for  suitable  optimization 
of  engine  bearing  failure  data. 

5.2  APPLICATION  TO  BEARING  FAULT  DATA 

Techniques  to  detect  and  analyze  the  approach  of  bearing  failure  in  rotating  machinery  have 
obvious  economic  benefits.  In  Fig.  2 1  are  shown  the  typical  Fourier  spectra  of  four  types  of  bearing 
failure.  The  spectra  are,  of  course,  obtained  after  failure  has  led  to  the  reestablishment  of  stationarity. 

The  signal  for  a  typical  impending  failure  is  shown  in  Fig.  22.  Note  the  long  interval  of 
stationarity  before  the  failure,  the  detailed  non-stationaiy  interval  at  failure  onset,  and  the  following 
stationary  interval  which  we  call  “stationary  failure.” 

Fourier  analysis  of  the  data  using  1024  point  transforms  is  shown  in  Fig.  23.  The  spectra  give 
accurate  information  for  the  two  stationary  regimes,  but  are  not  useful  in  the  transient  regimes. 
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The  Malvar  transform  of  the  signal  is  shown  in  Fig.  24  using  1024  point  transforms.  Even 
though  this  transform  is  cosine  based,  we  see  the  qualitative  similarity  for  the  three  regimes.  Indeed, 
spectra  for  the  stationary  regimes  ate  readily  recognized. 

Clearly,  a  long  Malvar  window  can  replace  the  shorter  windows  in  the  stationary  regimes.  In 
the  transient  regime,  on  the  other  hand,  much  narrower  windows  are  needed  to  display  significant 
details. 

Future  work  on  this  signal  will  include  the  implementation  of  the  basis  optimization  algorithm 
described  in  the  previous  section  so  that  the  splitting  and  merging  of  the  windows  can  be  done 
automatically  and,  at  some  point,  in  real  time. 

6.0  CONCLUSIONS 

Wavelet  theory  has  been  presented  and  contrasted  against  Fourier  analysis  from  a  basis 
approach  for  Malvar  and  Daubechies  wavelets.  The  wavelet  analysis  was  shown  to  perform 
similarly  to  that  of  Fourier  analysis  while  maintaining  the  element  of  time  in  the  transformation. 
The  theory  is  sdll  new  and  requites  additional  work  to  resolve  interpretation  issues,  but  it  has  been 
shown  to  be  a  valuable  analysis  tool  for  transient  data  analysis,  especially  if  combined  with  tradi¬ 
tional  Fourier  analysis  where  appropriate.  One  area  where  wavelets  offer  significant  analysis 
improvement  is  in  the  detection  of  transient  edges,  and  further  work  continues  in  this  area  for  image 
processing  and  other  applications.  New  wavelets  such  as  the  Malvar  and  Harmonic  wavelets  are 
becoming  available  which  have  promise  for  improving  the  application  and  interpretation  of  wavelet 
analysis. 

The  Malvar  wavelets  were  chosen  for  engine  bearing  failure  analysis  because  their  structure 
allows  representation  of  stationary  signals  typical  of  rotating  machinery  as  well  as  representation 
of  the  transient  effects  that  mark  the  onset  of  bearing  failure.  Their  similarity  to  the  Fourier 
transform  makes  them  familiar  to  engineers  already  accustomed  to  using  Fourier  spectra  to  under¬ 
stand  engine  behavior.  The  use  of  on-line  window  length  optimization  will  focus  attention  on  the 
unusual  behavior  of  an  engine,  while  more  normal  stationary  behavior  can  be  quickly  grasped  and 
assimilated.  Very  efficient  lossless  data  compression  is  an  immediate  benefit  of  the  optimization;  it 
also  generates  effective  pattern  recognition  parameters  for  use  by  higher  level  classifiers  such  as 
neural  networks. 

The  analysis  of  actual  engine  data  indicates  more  than  fair  likelihood  that  the  Malvar  approach 
will,  in  fact,  deliver  the  benefits  to  analysis  for  which  they  were  chosen.  All  but  the  optimization 
has  been  completed  with  algorithms  which  are  highly  efficient  and  readily  ported  to  digital  signal 
processing  computers;  these  algorithms  are  also  parallelizable  and  vectorizable.  Optimization  can 
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be  done  in,  at  worst  0(n  log(n))  time,  but  final  determination  has  not  yet  been  made  of  the  actual 
speed,  or  of  its  suitability  for  parallelization, 

Our  studies  of  both  pulse^echo  data  and  bearing  failure  detection  point  up  some  important 
facts  concerning  wavelet  analysis.  First,  wavelet  analysis  includes  a  wide  range  of  alternatives  and 
supplements  to  Fourier  analysis,  some  already  well  known  and  some  new.  As  such,  wavelets  do  not 
replace  Fourier  analysis;  they  broaden  it.  Second,  the  wavelet  technique  used  for  a  problem  must 
be  chosen  appropriate  to  the  problem  if  improvements  to  Fourier  analysis  are  to  be  realized. 
Furthermore,  where  Fourier  analysis  already  provides  well-understood  insight,  wavelet  techniques 
should  complement  the  Fourier  technique.  Third,  the  evaluation  and  use  of  wavelet  techniques 
requires  a  thorough  foundation  in  Fourier  analysis  as  well  as  additional  training  in  a  variety  of 
wavelet  analysis  techniques;  while  the  results  of  well-designed  wavelet  analysis  are  easily 
comprehended,  casual  acquaintance  with  the  techniques  is  simply  inadequate  for  the  judgment 
required  to  produce  a  successful  analysis  design. 
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Time,  sec 
a.  Example  signal 


Frequency,  Hz 


b.  FFT  spectrum  composite  signal  for  STFT  analysis 
Figure  1.  STFT  example  signal  and  FFT  spectrum. 
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a.  Windowed  signal  for  STFT  analysis  (FFT  Size  =  32) 


b.  Composite  spectrum  for  STFT  analysis  (FFT  Size  "  32) 
Figure  2.  STFT  spectrogram  example. 
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a.  Composite  signal  for  STFT  analysis  (FFT  Size  ■■  128) 


b.  Composite  spectrum  for  STFT  analysis  (FTT  Size  ■  128) 
Figure  3.  STFT  spectrogram  example  with  improved  high-frequency 
time  location. 
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NON-WINDOWED  Datafile.  Fig.  31 
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NON-WINDOWED  FFT  File:  Fig.  31 


WAVELET  Datafiles:  Fig.  33 
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Figure  6.  Comparison  of  Fourier  analysis  for  sine  waves  and  Haar  wavelet 
analysis  for  square  waves. 
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Figure  7.  DAUB4  wavelet  applied  to  STFT  example  signal  (Ref.  Fig.  14). 
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Figure  8.  DAUB  20  wavelet  applied  to  STFT  example  signal  (Ref.  Fig.  14). 
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a.  Wavelet  datafiles,  Fig.  59 


b.  Wavelet  datafiles.  Fig.  63 


Scale  (Frequency)  Index 


c.  Wavelet  datafiles.  Fig.  64 

Figure  9.  Wavelet  analysis  of  pulse-echo  simulated  data. 
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NON-WINDOWED  DataFile:  Fig.  59 


Sample  Time  Index,  int 


Frequency,  Hz 

Figure  10.  FFT  of  simulated  pulse-echo  data  of  Fig.  37a. 
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Sample  Time  Index,  int 


a.  Wavelet  detail  levels  10, 4,  and  3  only  for  Fig.  37a 
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b.  Wavelet  detail  levels  10, 4,  and  3  only  for  Fig.  37b 
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Sample  Time  Index,  int 

c.  Wavelet  detail  levels  10, 4,  and  3  only  for  Fig.  37c 
Figure  11.  Haar  wavelet  details  superimposed  on  simulatd  pulse-echo  signal. 
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WAVELET  Detail  Levels  10,  4,  and  3  ONLY 


WAVELET  DalaFiles  Fig  65 
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Figure  12.  Wavelet  analysis  of  pulse-echo  simulated  data  with  peak  not  coincident 
with  a  V  sample  point. 
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Overlapping  Malvar  Windows 
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Bell  Function  Profile  =  sin{it/4[l  +  sin(7i/2  ajt)]} 
where  P'{t)^+  |}’(t)"=  1 


ai'Kj  i^+Rj-1 

Figure  18.  Bell  ftinction  profile. 
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Figure  21.  Typical  Fourier  spectra  for  various  bearing  failures. 
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c.  Page  3 

Figure  22.  Time  domain  signal  for  impending  bearing  failure. 
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